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The use of APE smearing or other blocking techniques in lattice fermion actions can provide many 
advantages. There are many variants of these fat link actions in lattice QCD currently, such as FLIC 
fermions. The FLIC fermion formalism makes use of the APE blocking technique in combination 
with a projection of the blocked links back into the special unitary group. This reunitarisation 
is often performed using an iterative maximisation of a gauge invariant measure. This technique 
is not differentiable with respect to the gauge field and thus prevents the use of standard Hybrid 
Monte Carlo simulation algorithms. The use of an alternative projection technique circumvents 
this difficulty and allows the simulation of dynamical fat link fermions with standard HMC and 
its variants. The necessary equations of motion for FLIC fermions are derived, and some initial 
simulation results are presented. The technique is more general however, and is straightforwardly 
applicable to other smearing techniques or fat link actions. 
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Recent advances in computing power and in Lattice 
QCD (in particular, overlap fermions[l|) have allowed 
simulations at sufficiently light quark masses to see that 
the behaviour of quenched QCD can differ from the true 
theory both qualitatively and quantitatively in the chi- 
ral regime0,|3- As it is in the chiral regime where the 
difference from the quenched approximation will be high- 
lighted, we would like to simulate at light quark masses 
in dynamical QCD. This is an extremely computationally 
expensive endeavour. Ideally, this would be done using 
overlap fermions, although large scale dynamical over- 
lap simulations challenge the limits of current computing 
power, to say the least. 

Fat Link Irrelevant Clover (FLIC) fermions have shown 
a number of promising advantages over standard actions, 
including improved convergence properties [4] and 0(a) 
improved scaling without the need for nonperturbative 
tuning)^, [|| . Furthermore, a reduced exceptional config- 
uration problem has allowed efficient access to the light 
quark mass regime in the quenched approximation 
where recent studies have highlighted deviations from the 
true theory 13- As interest shifts ot focus on dynamical 
QCD, be it (truly) unquenched, or partially quenched, 
one might hope that the excellent behaviour at light 
quark mass displayed by FLIC fermions will carry over 
from the quenched theory to the unquenched one. This 
brings us to the issue of generating dynamical gauge field 
configurations with the fermionic determinant being that 
of the FLIC action. Brief accounts of this work were pre- 
sented last year 0, Recently, an alternative proposal 
for another type of smearing scheme that is differentiable 
has also appeared^. 



(HMC) [13. It is exact, ergodic and is 0(V*), that is, 
it scales almost linearly with the lattice volume V. In or- 
der to introduce our notation and a framework for our 
technique, we briefly review the HMC algorithm for gen- 
crating dynamical gauge field configurations. 

We wish to generate an ensemble {t/j} of (statistically 
independent) representative gauge fields distributed ac- 
cording to the probability distribution 



(1) 



where the effective action for full QCD 

S eS [U] = S g [U] - In det D { [U], (2) 
dct D f = [ ThpDi/) e ~ I d^i^D^ix) 



(4) 
(5) 



is obtained from the standard action 

S[U,$,tl>] = S s [U] + S f p,fail>], 

S| = / d*x i>{x)Di[U]vj{x) 



by integrating the fermionic degrees of freedom ip, ip of 
the functional integral 

(O) = -jr / VUV^Vyb 0[U,$,ip} e -s\P^M, (6) 



For Wilson-like fermions in the physical region (away 
from exceptional configurations) det Df is real and posi- 
tive. Hence if M = d\d u then det M = det Dj , and as 
M is positive definite the following integral is convergent 



I. HYBRID MONTE CARLO 
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detM" 



P^^-JA^WM-VW, (7) 



The standard technique for simulating dynamical 
fermions has for some time now been Hybrid Monte Carlo 



and as detM 



we have the fermion determi- 



dct A/- 1 

nant for two flavour QCD given in terms of an auxiliary 



2 



pseudo-fermion field 4>, so called because it is complex 
(bosonic) rather than Grassmannian. It is an identity 
that 



det Df = y dct D { det D\ = Vdct M = det VM, (8) 

and hence to simulate an odd number of sea quark 
flavours, it is possible to use det VM 11]. 

For HMC, the four-dimensional lattice theory is made 
five-dimensional through the introduction of a fictitious 
(simulation) time, or evolution parameter r. The gauge 
field U is associated with its (fictitious) conjugate mo- 
menta P, and the five-dimensional system is described 
by the Hamiltonian, 



H[U ) P]=Y d ^TrP M (.T 



S eS [U}. 



(9) 



with respect to the links |jy in order that the proposed 
configurations have high acceptance rates. Previously, it 
has not been clear how to perform HMC with fermion 
actions that make use of the APE blocking technique in 
combination with a projection of the blocked links back 
into the special unitary group. This reunitarisation is of- 
ten performed using an iterative maximisation of a gauge 
invariant measure, and this choice of reuniterisation is the 
source of the difficulty. The problem is that the iterative 
technique is not differentiable with respect to the gauge 
field and thus it is not possible to calculate Jj|, which is 
necessary for the equations of motion above. In the next 
section we consider an alternative technique and show 
that is does not suffer from this problem, allowing the 
simulation of dynamical fat link fermions with standard 
HMC (and its variants). 



For Gaussian distributed P the expectation value of an 
observable is unaffected by the 5D kinetic energy, 



(O) = ^ J VPVU 0[U]e- n ^ p \ 
Z = I VPVUe- n[u > p] . 



(10) 
(11) 



Given U, a new gauge field U' is generated by the update 
U — > U' , which consists of 

(i) Refreshment: Sample P from a Gaussian ensem- 
ble, p(P) cx e^^' TrP . Generate a pseudo-fermionic 
background field 4> according to p((f>) cx e~^ 1ju ^ . 

(ii) Molecular Dynamics Trajectory: Integrate Hamil- 
ton's equations of motion to deterministically 
evolve (U, P) along a phase space trajectory to 
(U>,P>). 

(iii) Metropolis step: Accept or reject the new config- 
uration (P',U r ) with probability p(U -> U') = 
min(l, e' AH ), AH = H[U'] - H[Uj. 

The discrctiscd equations of motion are derived by 
requiring that the Hamiltonian be conserved along the 
trajectory, = 0. The following discretised equations 
of motion then approximately conserve ri for small step 
sizes, At, 

U^x.t + At) = C/ Al (x,T)exp( i ArP M (x,r)), (12) 



P^x^ + At) = P„{x,t) - U„(x,t) 



5U M (x, t) 



(13) 



II. SU(3) PROJECTION 

The APE smeared links U^(x) present in the FLIC 
fermion action are constructed from U^{x) by performing 
n smearing sweeps, where in each sweep we first perform 
an APE blocking step, 



Ujt/J, 



w 



(14) 



followed by a projection back into SU(3),lljf\x) = 

P(y^\ x ))- Frequently, the projection is performed us- 
ing an algorithm which updates U^' iteratively in order 
to maximise the following gauge invariant measure, 

U^\x) S {U 1 S SU{3)\ReTr(U'V^ (x)) is maximal}. 

(15) 

We refer to this projection technique as MaxReTr projec- 
tion. While this projection minimises the local action|l4|. 
as we mentioned earlier it is not differentiable with re- 
spect to U/j,(x) and hence not suitable for use in HMC. 

Now, given any matrix X , then X'X is hermitian and 
may be diagonalised. Then it is possible (for det X =^ 0) 
to define a matrix 



W = X- 



1 



VxTx 



(16) 



whose spectrum lies on the complex unit circle and is 
hence unitary (w(z) = z/z*z is the complex version of 
the sign function). Furthermore, W possesses the same 
gauge transformation properties as X. This is easily seen. 
Let X^(x) transform as 



X^x)^G(x)X^x)GHx + e^, (17) 



In our implementation, we evaluate the matrix exponen- 
tial directly through diagonalisation, rather than expand- 
ing it. 

FLIC fermions are clover-improved fermions where 
the irrelevant operators are constructed using APE 
smeared links [l2l fl3| . As with other efficient updating 
algorithms, HMC makes use of the variation of the action Xj l (x)X fl (x) — > G(x + e^)XJ 1 {x)X ll (x)G\x + e M ). (19) 



then 



and hence 



XUx)^G(x + e M )XUx)&(x), (18) 
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Sweep Unit Circle 



MaxReTr 



0.866138301214314 

1 0.960313394813806 

2 0.980735000838119 

3 0.988384926461589 

4 0.992103013943516 

5 0.994182852413813 

6 0.995457365275018 

7 0.996293668622924 

8 0.996878305318083 



0.866138301214314 

0.960348747275940 
0.980751346847750 
0.988393707639555 
0.992107844842705 
0.994185532052157 
0.995458835653863 
0.996294454083006 
0.996878710433084 



TABLE I: The mean link uo for a single configuration as a 
function of number of APE smearing sweeps at a = 0.7, for 
the two different projection methods. The boldface indicates 
significant digits which match. The configuration is a dynam- 
ical gauge field with DBW2 glue and FLIC sea fermions, at 
(3 = 8.0, k = 0.1280 (a = 0.17 fm, m q = 200 MeV). 



Noting that [Xj 1 (x)X^(x)] 2 has the same transforma- 
tion properties as Xj l (x)X^(x) it is then straight forward 
to see that 



W^x) -> G{x)W^{x)G\x + e„), 



(20) 



as required. 

Given the unitary matrix W, we can construct another 
matrix, 



W 



1 



VdetW 



-.W 



(21) 



which is special unitary. Earlier work|15j has incorrectly 
omitted the cube root. As there are three different com- 
plex roots, we have a Z3 ambiguity which we break by 
choosing the principal value of the cube root 1 In select- 
ing the principal value, the projected matrices lie closest 
to those given by the MaxReTr method, and are hence 
smoother. The mean plaquette is closer to unity thus 
minimising the action. We refer to this technique for 
projecting X^(x) into the special unitary group as unit 
circle projection. 

The two methods produce smeared links that are differ- 
ent but numerically close (according to the usual matrix 
norm, = \J \ max {A^A)). Using the mean link as a 
measure of the smoothness of the smeared gauge field, 
Table [I] indicates that the two methods presented here 
produce equally smooth gauge fields. 

While numerically the two methods may be nearly 
equivalent, unit circle projection possesses a significant 
advantage over MaxReTr projection. The matrix in- 
verse square root function can be approximated by a 



1 For complex z, the principal value of the cube root satisfies — ^ 
arg ^fz < ^. For purely real 2, we choose ^fz to be real. 



rational polynomial (whose poles lie on the imaginary 
axis) [H 03, W[X] « W k [X], 



W k [X} = d X(X^X + c 2n )J2 



1=1 



X^X 



(22) 



where the formula for the coefficients cio,6;,C; can be 
found Ref. |17| . This approximation is differentiable in a 
matrix sense for all X for which the inverse square root 
can be defined. This means that we can construct |jy for 
fermion actions which involve unit circle projection, and 
hence it is a reuniterisation method which is compatible 
with HMC. 



III. EQUATIONS OF MOTION 

Having now defined the APE smearing prescription 
(with projection) in a differentiable closed form, we pro- 
ceed to derive the equations of motion necessary for the 
use of the HMC algorithm with FLIC fermions. 



A. Mathematical Preliminaries 

The equations of motion are derived using multi- 
variate calculus. To make the derivation simple and 
provide an understanding of how best to implement the 
equations efficiently, we develop some appropriate math- 
ematical tools. Using index notation, we define a (min- 
imal) set of tensor operations (including differentiation) 
such that we can perform the derivation in an index free 
language. 

The derivative of a real- valued function f[A] with re- 
spect to the matrix A is a rank 2 type (1,1) tensor (dis- 
tinguishing contravariant and covariant indices), 



r - 
[ dA l3 



d 



■f[A]. 



(23) 



The derivative of a matrix- valued function M [A] with 
respect to the matrix A is a rank 4 type (2,2) tensor, 



dM 
[ ~dA^ 



d 



3 1 



■M[A] 



(24) 



The set of type (m, n) tensors T™ forms a vector space. 
We define the outer product <g) : T\ x T\ — > T\ 



as 



(A ® B) 1 j k i 



A 1 3 B 



(25) 



Noting carefully the index ordering, define the "direct" 
product © : T\ x T\ -> T\ as 



{A@B) lk 



j 1 



A jB l 



(26) 



Given a scalar function f[B] and a matrix function B[A] 
the (scalar-matrix) chain rule states 



dA 



8f_ dB_ 

dB*dA : 



(27) 
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where we define the contraction induced by the chain rule The (matrix-matrix) product rule is 
as the (rank 2) star product, * : T\ x T\ — > T\, with 



(A-kTYi = A 3 k T' 



i i- 



(28) 



Given two matrix functions M[_B] and -B[A], the (matrix- 
matrix) chain rule states 



dM 
~dA 



dM dB 
~8B *0A' 



(29) 



where we define the contraction induced by this chain 
rule as the (rank 4) star product, * : T§ X T§ — ► T|, 



with 



(S*Ty j k l = S t m n l T 



m k 
j n- 



(30) 



It is interesting to note that the star product induces an 
algebra structure on the vector space of type (2,2) ten- 
sors, that is, (7~2,+,*) is an algebra with multiplicative 
identity I (g) I. 

We define juxtaposition for A e T\,T e l\ by the 
contractions 



(AT) 1 k 

(TAy 1 



3 I 
k 
3 I 



T 



i k 



3 U 

A m t. 



(31) 
(32) 



All our derivatives will be derived from the basic ma- 
trix differentiation rule. Given matrices M, A, B, C then 
for M = ABC we have 



dM 
~dB 



A®C. 



An immediate consequence of this is that 
dM 



dM 



= Ji 



The (scalar-matrix) product rule is 



(33) 



(34) 



(35) 



d ( XY )-X dY + dX Y 
dA {XY) - X dA + dA Y ' 

which is easily shown to imply the identity 



dX- 



dA 



(36) 



(37) 



In the following sections we will make use of the identity 



A + {B®C) = {A-B)C, 



(38) 



where A- B — A l jB J i. Additionally, of particular numer- 
ical importance is the identity 



A * (B <g> C) = BAC, 



(39) 



which has two major benefits. It allows us to evaluate 
two matrix multiplications instead of an outer product 
(computational saving), hence enabling us to implement 
the equations of motion without having to store any ten- 
sor fields (memory saving). 

B. Standard Derivatives 

The equations of motion for FLIC fermions are de- 
rived starting from the equations for the standard clover 
fermion action. We divide the effective action into its 
gauge part and pseudo-fermionic part, 



S e g = Sg + 5, 



(40) 



We reformulate some standard results in terms of the 
mathematics of the previous section. We will adopt a 
more convenient notation for quantities with a lattice site 
index x, using a subscript U^x rather than U^x). The 
matrix products of link variables are often denoted dia- 
grammatically. 



For a plaquette plus rectangle improved gauge action, 

S g = J2 ReTr[/3 lxl (l-C/ AlI ,(.T)) + /3 2xl (l-i?^ 1 (.x)) + /3 lx2 (l- J Rif (x))] , 



(41) 



we have 



dSe 



dU, 



ft ■ X 



S A - ( n+ u) 



+ ^2x: 



(42) 



where the filled circles indicate the point x. The coefficients /3 n 



\f3c mxn depend on the choice of gauge action. 



For Liischer-Weisz glue, ci x i = §,C2 X i = ci X 2 
such that cixi = 1- 

The pseudo-fermionic action is S p { 



j2- For DBW2 glue, we choose coefficients that are normalised 



J2x 't'ivx, where r\ = (D^D) 1 <fi, hence by equations ll'UiH and (|37|l we have 



■5 



9S 2L = ^tpt^-i ( D l™ + 1EL D ) (D t D) - V . (43) 



Setting x = -D77, we obtain 



7? + ^— — x. (44) 



dU^ x dU^ x dU^x 
Now, the FLIC action is explicitly given by 

1 C/ fl U U R1 > W 1 

(^flicV)x = - 7^)^ + (-Jr + + (4 + m - a^^JV-x (45) 

and contains three terms, the Dirac term (constructed with standard links), the Wilson term and the clover term 
(using fat links for Fj* x , and setting = 5 [7^, 7^] ) ■ Hence we may decompose the pseudofermionic derivative into 
three terms also. The first comes from the Dirac term, 



dS^ _ 1 



Tr spin (r7| <g> l^Xx+n + xl® 7nVx+»), (46) 



dU^x 2u 

while the Wilson and clover terms only explicit dependence is on the smeared links, 



dS w i 

ds;\ 



where the vector outer product defines a matrix (77 (g> = ^'Xj- The one loop clover term is given by F v \ <y = 

l{C v \, y - Cl X y ), where 

C v \, y = U( +u ) ( +A)i y + U(+X)(- V ),y + U(- X )(+v),y + U (-v)(-\),y, (49) 

and Ut± v -\i±\\ y indicates the plaquette starting at y, oriented in the v — A plane, with the first (second) link in the 
direction indicated by the first (second) index. When differentiating with respect to x , any terms where y lies 
further away from x than x ± v ± A will be zero. Further, noting that the derivative is zero unless either v = fi or 
A = /i and v ^ A we can without loss of generality choose /i — v. Letting fi — v be in the horizontal direction and A 
be in the transverse (vertical) direction, the contribution to the derivative due to the clover term is 



dF 

' 'L-J-. — _ I / 1 I V'/< ... J I A - / A ' ." _ 

J 

where the filled circles indicate the point x, and the point y is located at the start (end) of the diagrams that lie on 
the left (right) side of the outer product, as can be deduced from the Kronecker S-s. 



Smeared Link Derivatives 



Now, having constructed the explicit derivatives of S p f 
with respect to the thin and fat links, the total derivative 



with respect to the thin links is 

dS pf _ dS p( 8S pi 



dU?, 



dU u 



dU u 



dU u 



(51) 



If we have performed n sweeps of APE smearing to form 
the fat links, then the right hand term is constructed 
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through n applications of the chain rule 



Using 



dS OS dU$ , dS duj$ 



_ _ i _ * /g2) 



until we arrive at 



dU 



We note here that the partial 

derivative with respect to a forward (backward) link only 
picks up terms that contain the forward (backward) link, 
and not its conjugate (that is, U and W are consid- 
ered independent with regard to partial differentiation, 
see Eqs. l(T3 )l .(|5S |l .l(57 )l l. For the sake of both compu- 
tational efficiency and simplicity, this chain rule is itself 
composed of several chain rules, and hence evaluated in 
several steps. Each step corresponds to a step in the APE 
smearing process, but we go through them in reverse or- 
der. 

The final step in the APE smearing process (with unit 
circle projection) is 



[/(") = w {n) 

detWXaN 



(53) 



Therefore the first chain rule corresponds to this step, 



dS dU, 



(n) 



dS 



fl.X 



detW.J/®/ J, 



where 



det A = e ijk A 1 i A 2 j Al 



(54) 



(55) 



and hence denoting the permutations 7T; = (i mod 3) 

1, 7T 2 = TIV. i 



d det A 



— e ; A W < , A^i 



(56) 



There are several chain rules that correspond to 

W(n) = y(n)(y(n)1y("h-i (57) 



For the first, we define H^ ^ = V/lx^n^- Then 



as aw u 



OS 



dH^x dW^x dH^ : 



(58) 



W^.x « d Q V^ x {H^x + c ) ^2 



(59) 



we have 



dW^x , T . ( T<> xr b l 



dH, 



jt ..c 



1=1 



Hu,x + ci 



k 

(H^x+co^k-— 



i=i 

We can then construct 



Hn.x + Q Hn 



ci 



(60) 



as ^ dW ux dS 8H UX 

— * — — h ~ * — 



dV^x dW^x dV^x dH^x dV^ 



9S IT „-i. dS , TA+ 

*(I®H tt ») + -^ r -*(vL x ®I), 



dW, 



(61) 



and also 
dS 



dS awt , 95" 8H„ X 



dVl x 



dWlx dVlx dH^x OVlx 

dS -i 55 

* ® /) + TTTj— * (/ ® V^, x ). 

f-L,X 



Last, we make use of the chain rule 
where 

V^ !3/ = (l-a)^+|^n + 

v^,, = (i - + , . , 
6 u 

and y is indicated by the filled circle. 



(62) 



x) i (63) 



(64) 
(65) 



It is then straightforward to show that 



dU, 



(n-l) 



= (1 - a)/ 



6^ 



I KUx-u- (66) 
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and 

dV,t, a 



where in these diagrams the filled circle indicates the point y. Hence, 
dS , . dS dS f I dS 



— "--E 1 ®r + T ® 1 C+A> (67) 



■7-— — ~k I <Si j + — — * . <g>Z + 1 + 1 * . ®J, 68 

9F„ aK— ■■ ■■ T ^ T/t ^ T/t 1 



where now x is indicated by the filled circle. 

Having constructed the total derivative of the action with respect to f7 ftx , we can calculate the variation of S with 
respect to the gauge field (noting W — U" 1 ), 

su - dU + dW du - dU dW ' 1 j 

and hence the necessary equations of motion, (|12|l and (|13|) . It is numerically efficient to make use of the fact that 



dS /dS xt 



dW V dU 



(70) 



IV. SIMULATION RESULTS 

We have implemented the equations above within a 
standard two-flavour HMC, with multiple time scales. 
Expensive pseudo-fermion momenta updates are per- 
formed at a larger step size At = Ar p f and the 
cheaper gauge momenta updates are performed more of- 
ten, Ar g = iAr p f, for some integer n. Molecular dy- 
namics trajectories are of unit length, n m( jAr = 1. In 
particular, we have implemented a modified version of 
the Ritz algorithm to diagonalise arrays of 3 x 3 matrices 
in parallel. This routine is used in the SU(3) projection 
step, and is also used to calculate the matrix exponentials 
that are needed in other parts of the algorithm, avoiding 
the need to use polynomial approximations to the expo- 
nential. This means that the accuracy of the exponential 
in Eq. l(T2|) does not depend upon the step size At. 

An eighth order Zolotarev approximation to the in- 
verse square root is used to approximate W^x in unit 
circle projection. We find that the spectral range at this 
order is ample. In smooth gauge backgrounds it is easily 
shown that unit circle projection is well defined, that is 
detV^V^ > 0. 

If we assume a smoothness condition ||1 — L^i^H < 
e Vx, /i, is, then we have a lower bound 

Vl x V^ x >l-2ae- a 2 e 2 . 

To prove this, we note that APE blocking may be written 



in terms of the plaquette field, 

= U^ x (1-2^(1- uU) . (71) 

Define Z = £ ±I ^(1 - U^ x ), then 

Vl x V„ x = 1 - ~(z + Zt) + £z*Z. (72) 
As \\Z\ \ < 6e, we then have 

Amin^ > 1 - 2ae - aV, (73) 

which is strictly positive for small enough e. 

While the smeared link equations of motions are com- 
plex, our implementation evaluates them efficiently due 
to the optimisations that can be performed through 
the calculus we constructed earlier. At large sea quark 
masses the code already spends over 90% of its time in the 
CG inversion required to calculate 77 — (D^ D)^ 1 ^, and 
as the quark mass decreases this fraction increases. So as 
is standard, the generation of dynamical gauge fields is 
dominated by the CG inversion. Simulation results are 
presented in Table ITU 

It is a simple exercise to apply our results to gener- 
ate gauge fields with dynamical FLIC Overlap quarks, 
although this would be extremely computationally inten- 
sive. The availability of HMC as a simulation algorithm 
for dynamical FLIC fermions is significant, as it scales 
almost linearly with the lattice volume V, whereas pre- 
viously there were only 0(V 2 ) alternatives [lsj. Fur- 
thermore, the method we have described is general and 
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can be straightforwardly applied to any fermion action 
with reuniterisation, including overlap fermions with a 
fat link kernel [3 S El E3 , or other types of fat link 
actions l23l that may involve alternative smearing tech- 
niques 24] . Additionally, any of the variants of HMC 
can be also be used, in particular Polynomial HMC|25| 
or Rational HMC|2(| which allow for the simulation of 
odd numbers of sea quark flavours. 
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TABLE II: Simulation parameters and results for various dynamical simulations. The parameters are the gauge coupling, 
hopping parameter, gauge action, step size, and psuedofermion to gauge step size ratio. The results given are the mean link, 
lattice spacing (in fm, obtained from ro via the static quark potential) and pion mass (in GeV). Two degenerate flavours 
of FLIC sea quarks are used, with either Lusher- Weisz (IMP) glue or DBW2 glue. These results are obtained from 20 
12 s x 24 configurations. Simulations are done using multiple time step HMC with trajectories of unit length. Although an 
exact comparison is difficult, for a given step size and quark mass the acceptance rates obtained compare well with standard 
simulations (see for example Ref. |27fl). 



